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Abstract 

We report on a method for matching the next-to-leading order calculation of QCD jet production 
in e + e~ annihilation with a Monte Carlo parton shower event generator (MC) to produce realistic 
final states. The final result is accurate to next-to-leading order (NLO) for infrared-safe one-scale 
quantities, such as the Durham 3-jet fraction y%, and agrees well with parton shower results for 
multi-scale quantities, such as the jet mass distribution in 3-jet events. For our numerical results, 
the NLO calculation is matched to the event generator Pythia, though the method is more general. 
We compare one scale and multi-scale quantities from pure NLO, pure MC, and matched NLO-MC 
calculations. 
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I. INTRODUCTION 



Perturbation theory is the basic tool for deriving predictions for experiment from the 
standard model and its extensions. When some of the particles feel the strong interaction, 
a perturbative expansion in powers of the strong coupling a s can be employed. As long as 
a short-distance process is involved, the reasonably small value (a s ~ 1/10) of the running 
strong coupling justifies such an expansion, and the leading term is a good first approxima- 
tion. However, the estimated error to a leading order (LO) prediction (that is, the error that 
one estimates from leaving out all of the higher order terms) is often so large that one wants a 
next-to-leading order (NLO) calculation. In fact, next-to-leading order predictions are avail- 
able in the form of computer programs for a wide variety of processes for electron-positron 
collisions, electron-hadron collisions, and hadron-hadron collisions. The essential limitation 
is on the number of partons involved in the hard interaction. Thus p + p — > W + 2 jets + X 
is available at next-to-leading order while p + p — > W + 3 jets + X is not. 

A typical next-to-leading order calculation predicts the expectation value (S) of an 
"infrared-safe" observable S in the style of an event generator. That is, the program gen- 
erates a large number N of events characterized by their final states /„, with each event 
having a weight w n . Letting the value of S for state /„ be S(f n ), the calculated expectation 
value of S is 

(s) = ^j:w n s(f n ). a) 

JV n=l 

(A convenient normalization is (i/N)^2^ =1 w n = 1 so that (1) = 1. One could also take 
(1) = (Jtot-) The quantity (S) has the perturbative expansion 

(S) = C (S) of + of +1 + C 2 (S) of + 2 + • • • , (2) 

where B is the power of a s in the lowest-order graph for (S) . The numerical result of the 
NLO calculation is not exact, but produces the first two terms of this expansion. 

Typical programs that do next-to-leading order calculations suffer from a serious defi- 
ciency: the final states /„ are not realistic. First of all, they consist of just a few partons. 
In the case of e + + e~ ^3 jets, which we consider in this paper, the final states consist of a 
quark, an antiquark, and one gluon or else a quark, an antiquark, and two gluons or another 
quark- ant iquark pair. This is not the kind of state that can be used directly as input to a 
detector simulation. 

The situation is even worse than might be deduced from the fact that the final states 
consist of few partons instead of (many) hadrons. At its core, an NLO calculation depends 
on the use of an infrared-safe inclusive observable since it depends on a cancellation of 
divergences associated with different numbers of on-shell partons. As a concrete example, 
consider the fraction of events in e + + e~ — > hadrons that contain precisely three jets as 
determined by the kx (or Durham) algorithm [1] with y cut = 0.05. Call this fraction f 3 . We 
have 

(/ 3 ) = C (f s ) a s + d(/s) a 2 s + C 2 (/ 3 ) a 3 s + ■ ■ ■ . (3) 

A next-to-leading order program can accurately predict C and C\. Now suppose that 
we use the program to determine the differential three-jet fraction dfe/dM, where M is 
the mass of one of the jets (so that each event contributes three entries for dfs/dM and 
JdM dfe/dM = 3/3). Even if one is ultimately interested in / 3 , it is reasonable to be in- 
terested in predicting this more detailed quantity, because a realistic detector may respond 



2 




5 10 15 20 

M [GeV] 



FIG. 1: Jet mass distribution in three-jet events, f$ dfe/dM, calculated at next-to-leading order 
for y^s = Mz- The three jets are identified using the algorithm with y cut = 0.05. Then M is the 
invariant mass of one of the three jets, with each event making three contributions to dfe/dM. The 
renormalization scale is chosen as fj, = and a s {Mz) = 0.118. This is a pure NLO calculation 

using [5]. There is a large negative contribution in the first bin. 

differently to narrow and wide jets and we may be concerned with how much detector effects 
influence the measurement. 1 Now, one can expect the NLO result to be approximately cor- 
rect for dfs/dM as long as M is large, but for Mj^fs < 1, a sophisticated user will expect a 
pure NLO program to give an answer that is not physically realistic. We illustrate what goes 
wrong by simply plotting the result for f^ 1 dfe/dM versus M from a purely NLO program 
in Fig. 1. We see that f^ 1 dfe/dM becomes large for small M and that the fraction of events 
in the first bin is large and negative. As M —>■ 0, df 3 /dM according to the NLO calculation 
develops a log(M)/M singularity together with a term proportional to 5(M) with an infinite 
negative coefficient. The integral /3 = (1/3) JdM dfo/dM of this highly singular function 
gives an accurate estimate of the three-jet fraction f 3 but the differential distribution is 
completely unphysical. In fact, it requires a binning of the distribution even to get finite 
numbers. A much more realistic expectation for f^ 1 dfe/dM can be obtained by using a 
parton shower Monte Carlo event generator [2-4]. Recall that the shower Monte Carlo pro- 
grams add up a series of logarithmic terms, approximating the behavior of Feynman graphs 
near the singularities corresponding to collinear parton splitting or soft gluon emission. The 
prediction of the shower Monte Carlo Pythia [2] for f^ 1 df^/dM is shown in Fig. 2 and is 
compared to the NLO calculation. We see that the NLO perturbative calculation agrees 



Actually, for this purpose one would be more interested in dfe/dAIi dM-i dM^, but we discuss df^/dM for 
reasons of simplicity. 
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FIG. 2: Jet mass distribution in three-jet events, f% dfs/dM, calculated according to Pythia [2] 
with default parameters for = Mz- The cross section is defined as in Fig. 1. The jet mass 
distribution from the perturbative NLO calculation in Fig. 1 is shown for comparison. 

qualitatively with Pythia for large M, but fails for small M. 

Clearly if one were interested in how the differential response of a detector to wide 
and narrow jets would affect a measurement of f%, one would be better off using Pythia 
than the NLO program, even though Pythia has only leading-order accuracy for f% = 
(1/3) / dM dfz/dM. However, would it not be better to add parton showering and hadroniza- 
tion according to Pythia to the NLO calculation? That is the topic discussed in this paper. 

There has been substantial recent progress in adding a parton shower and hadronization 
to an NLO calculation. Frixione, Nason and Webber [6-9] have presented and implemented 
an algorithm that makes use of the shower Monte Carlo program Herwig [3] to do NLO 
calculations for a selection of processes for which there are two hadrons in the initial state 
and at leading order the final state particles that carry color are massive. Processes available 
include the hadropro duct ion of single vector and Higgs bosons, vector boson pairs, heavy 
quark pairs, and lepton pairs. Two of the present authors have described an algorithm for 
a process with massless color carrying particles in the final state, namely e + + e~ — > 3 jets 
[10, 11]. This algorithm was implemented with an NLO calculation coupled to a small 
self-contained shower generator with no hadronization. 

The calculations of Refs. [10, 11] sufficed to demonstrate the principles of the algorithm, 
but for practical purposes one would like to have a shower algorithm that has been tested 
against data (as in Herwig and Pythia) and one would like to include hadronization. For 
this reason, we present here results from coupling the program described in Refs. [10, 11] 
to the standard shower Monte Carlo program Pythia, including its hadronization following 
the Lund string model. The code described here is available at the site [12]. Other work on 
this general subject can be found in Refs. [13]. 
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FIG. 3: Parton splitting in Refs. [10, 11]. The filled circles represent graphs for the Born amplitude 
and complex conjugate amplitude. Each of the partons emerging from the Born amplitude splits 
into two partons with a vertex, represented by the squares, that includes a Sudakov suppression 
factor. The extra gluon coming from the circles on the right and left represents the soft gluon 
radiated from the three jets. Each of the seven daughter partons undergoes further, secondary, 
splittings according to Pythia and enters the final state as a complete shower with hadronization. 
The secondary splittings are represented by the diamonds. 

II. SKETCH OF THE ALGORITHM 

The general algorithm that we use to couple an NLO program for e + + e~ — > 3 jets to 
a shower Monte Carlo is described in some detail in Refs. [10, 11], so we present just a 
sketch here. Specific details about the coupling to Pythia that go beyond Refs. [10, 11] are 
described in the following section. 

Consider first the Born graphs. A quark-antiquark-gluon state is generated with a weight 
proportional to one of the Born graphs times the complex conjugate of one of the Born 
graphs. Each of these partons splits into two with a probability described by a splitting 
function and a Sudakov exponential that represents the probability not to have split at 
a higher virtuality. This splitting is called a primary splitting and is represented by the 
square vertices in Fig. 3. The splitting functions have the right singularities to represent 
the collinear limit of QCD matrix elements and the collinear x soft limit, but they are not 
correct in the limit of the emission of a wide angle soft gluon. Accordingly, we generate 
a gluon emission with a weight that matches the matrix element for the radiation of a 
very soft gluon from the antenna produced by the three outgoing partons. This emission is 
represented by the gluon in Fig. 3, which is drawn to suggest that it is emitted from the 
outgoing partons as a whole but not from any particular one of them. 

Next, the seven parton final state thus generated is modified somewhat as described in 
the following section and fed to Pythia, which creates a complete shower and hadronization. 
We have modified Pythia slightly so that it is able to accept a partially developed shower 
and generate suitably narrow jets from each of the seven partons. We outline the required 
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FIG. 4: Treatment of order af + graphs. One particular cut diagram is shown, in this case a cut 
diagram with a virtual loop. The second diagram illustrates a subtraction term, which corresponds 
to a contribution at order af +1 from Fig. 3. The partonic final state is generated with a weight 
proportional to the difference of matrix elements, then this final state is sent to Pythia, which is 
represented by the diamonds. 

modifications in the following section. 

The splitting incorporated in the Born graphs has two important features. First, when the 
virtuality q 2 of the pair of daughter partons is small, the splitting probability is proportional 
to (P(x)/q 2 )dq 2 dx where x is the share of the momentum carried by one of the daughter 
partons and P(x) is the appropriate Altarelli-Parisi splitting function. Second, the collinear 
singularity at q 2 — > is damped by a Sudakov factor with the behavior exp(— a s clog 2 (g 2 )) 
for q 2 — > 0. Thus each of the three original partons makes a jet, but in the limit of small 
a s , each jet is usually very narrow and appears in an infrared-safe measurement like a single 
massless parton. There is a qualifying adverb "usually" here. A fraction a s of the time, one 
of the splittings has a substantial virtuality and we get a four-jet final state. Thus there is 
an order af +1 effect in which some probability is removed from the three-jet final state and 
given to a four-jet final state. In a purely NLO calculation, this effect is included in the af +1 
graphs. Since we do not want any part of the order af +1 contribution to be counted twice, 
we need to subtract these terms from the af +1 graphs. One of these subtraction terms is 
illustrated in Fig. 4. The subtraction terms have the effect of removing collinear and soft 
divergences from the order af +1 graphs. 

In addition to the Born graphs, there are subtracted order ctf +1 graph with either three 
or four partons in the final state. Each of the final state partons from these graphs is fed to 
Pythia, which creates a complete shower and hadronization using appropriate initial condi- 
tions as described in the following section. This showering is represented by the diamonds 
in Fig. 4. 

The Sudakov suppression factors play an important role in the primary showering that 
takes place before the partons are passed to Pythia. In the case of quark splitting, the 
Sudakov factor has the form 

where q 2 is the virtuality of the splitting and q is the momentum of the parent quark. The 
function V g / q (J 2 , z) represents "virtual" splittings and is derived from the one loop quark self- 
energy graph in the Coulomb gauge. It reduces to the usual Altarelli-Parisi splitting function 
in the limit I 2 — > at fixed z. However the usual 1/z singularity is absent, being cut off for 
z < l 2 /\q 2 \, where q is the momentum of the mother quark. The integration over z produces 
a logarithm of l 2 /\q 2 \. Then the integration over I 2 produces a second logarithm so that 
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splitting at very small q 2 is suppressed by a Sudakov factor that is, in a first approximation, 
exp(— a s clog 2 (g 2 )). In Refs. [10, 11], the Sudakov factor is not derived as a summation of 
logarithms. Rather, it is derived as a useful way of expressing the effect of virtual graphs 
on small q 2 splittings. Compare this to the standard perturbative treatment of small q 2 
splittings in which large positive probabilities from real parton emissions are counterbalanced 
by large negative probabilities from virtual parton graphs. With the Sudakov factor, we avoid 
the small q 2 singularities: we suppress them with an exponential factor. When expanded 
in powers of a s , the result is the same as given by the standard perturbative calculation. 
Next, one simply notes that the Sudakov suppression factor is just what one wants as the 
first step in a shower Monte Carlo algorithm, where this factor appears as a way of summing 
logarithms. 

This algorithm is based on the order af and af +1 graphs for e + e _ — > 3 jets. It is designed 
to produce NLO calculations for infrared-safe three-jet observables while, at the same time, 
providing a reasonably accurate parton shower description for the inner structure of the three 
jets. Future NLO-MC hybrid programs will likely have a mechanism for switching from a 
three-jet description to a two-jet description, as in Refs. [14] and [15] at lowest order and 
Ref. [16] at NLO. Such a mechanism is not built into the present program. For this reason, 
the present program should not be used to calculate observables that receive significant 
contributions from two-jet final states. In fact, in the program, events that are too close to 
a two-jet configuration based on the thrust of the perturbative event are simply cut from 
the calculation. 

III. MERGING WITH PYTHIA 

In this section, we address the issue of how to pass the perturbative events after pri- 
mary showering through a standard event generator to include additional jet structure and 
hadronization. In this paper, we use the Pythia event generator, though there is no limita- 
tion (in principle) to using another one. 

For some time now, Pythia has been able to perform parton showering and hadronization 
on externally-generated partonic configurations. This has become technically simpler with 
the adoption of a general data structure by several Monte Carlo authors [17]. The parton 
shower that is added, however, is based on partons that are generated from a hard process 
at a single scale. We need something different since the further showering of the partonic 
state sent to Pythia needs to recognize several scales that were involved in generating these 
partons. 

We begin with seven partons from a Born graph as depicted in Fig. 3, three partons 
from an order af +1 graph with a virtual loop as depicted in Fig. 4, or four partons from 
an order af +1 graph with no virtual loop. These final state partons are generated with 
flavor identifications {g, u, u, d, . . .} and with labels indicating color connections based on 
the color structure of the underlying Feynman graphs. (The idea of color connections as 
used in parton shower Monte Carlo programs neglects terms suppressed by l/N 2 , where N c 
is the number of colors, and neglects quantum interference graphs. Thus the assignment 
of color connections can only be a rough approximation in a calculation with exact color 
factors for the graphs and with quantum interference included. Nevertheless, an assignment 
of color connections is needed in order for Pythia to generate string hadronization.) We 
also supply the quark masses used in Pythia to the final state quarks. The momenta for 
final state partons were generated in the approximation that all partons are massless, but 
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the momenta are adjusted so that final state quarks are on-shell with the proper masses. 

Before the partonic final state is passed to Pythia, a certain amount of analysis is needed. 

The first step is to account for the minimum virtuality of 1 GeV 2 allowed in a Pythia 
shower. In contrast, there are no infrared cutoffs in the virtualities of the parton splittings 
in Fig. 3 or in the energy of the soft gluon there. Thus in order to generate partonic states 
that are sensible within the context of Pythia, the program checks whether any of the 
three splittings indicated in Fig. 3 had a virtuality less than the minimum. If there was an 
unallowed splitting, the two daughter partons are replaced by the mother parton and the 
parton is flagged so that it is not showered by Pythia. In addition, if the soft gluon in Fig. 3 
(or any other final state gluon that did not arise from a 1 — > 2 splitting) has an energy less 
than 100 MeV, the gluon is erased and its momentum and color is reassigned to the other 
partons. This step leaves us with seven or fewer partons to be passed to Pythia. 

The next step is to generate a synthetic shower history that could have resulted in the 
final state partons. This shower history is generated based on the /c^-algorithm modified to 
account for the flavors and colors of the final state partons, so that we get a shower history 
that is allowed in QCD at leading order in 1/N%. In generating this shower history, we define 
the resolution variable to be the LUCLUS resolution variable that is used in Pythia: 

^ = 2(1 -003^)^^. (4) 

Here the energies and angles are defined in the overall center-of-mass frame of the event. 
The jet algorithm examines all of the kT,ij values for pairs {i, j} of partons that have flavors 
and color connections that would allow them to be combined. The pair with the smallest 
kr,ij is combined (by adding their four- momenta). Then the algorithm repeats until, at 
last, we have just one remaining quark-antiquark pair. We use this synthetic shower history 
to define, for each final state parton i, a certain initial virtuality scale Qi, a maximum 
transverse momentum A;™ ax , and a maximum angle #™ ax . The first splitting generated by 
Pythia for parton i should have kr < &™ ax and 9 < #- nax - 

To define /c™ ax , we find the kr of the splitting that produced parton i in the synthetic 
shower history. We define /c™ ax to be the lesser of this and a global maximum A>r, which is 
the largest kr,ij in the shower. The purpose of using this global scale is to guarantee that 
the hardest interaction is generated in the perturbative NLO calculation, not by Pythia. 

To define 6 l ° iax , we find the one or two final state partons j with which parton i is color 
connected and define #™ ax = min,, 9^. 

To define the initial virtuality scale Qi there are a number of cases. When parton % is a 
quark or antiquark that can be traced back to the 7* jZ vertex, we define Qi = s. When 
parton i is a quark or antiquark that can be traced back to a g — * qq vertex, we define 
Qi = Pg where p g is the four-momentum of the gluon at which the quark or antiquark line 
started. When parton i is a gluon that originated at a q — > qg or q — > qg vertex, we define 
Qi = Pra where p m is the four-momentum of the quark or antiquark that is the mother for 
the gluon. When parton i is a gluon that originated at a g — > gg vertex and is the less 
energetic of the two sister gluons, we define Qi = p^ where p m is the four-momentum of the 
mother gluon. When parton i is a gluon that originated at a g — > gg vertex and is the more 
energetic of the two sister gluons, we define Qi = p mm where p mm is the four-momentum of 
the grandmother of the gluon. 

The idea now is to ask Pythia to generate further showering for each final state parton 
i according to the standard Pythia splitting kernels, but with the first splitting limited by 
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the requirements kr < /c™ ax and 6 < 0™ ax . In order to incorporate the kx restriction into a 
shower based on the initial virtuality scale Qi, we follow the procedure of [15], based on the 
suggestion of [14] . The shower for each parton i is started at scale Qi but splittings that do 
not satisfy k T < k™ x are ignored ("vetoed"), where k T is the parton shower definition of 
this quantity, k\ = z{l — z)m 2 , using z as the energy fraction of one daughter with respect 
to the mother of virtuality m. Splittings that do not satisfy < #- nax are also vetoed, as 
is standard in Pythia. Once a first allowed splitting has been generated, the rest of the 
shower from parton i is generated according to the normal algorithms of Pythia. After 
parton showering, Pythia generates hadrons by string fragmentation. 



IV. RESULTS 

In this section, we test the program that combines the NLO calculation with showers and 
hadronization as sketched in Sees. I and II. We shall refer to this program as NLO+PS+Had. If 
we use just the Born graphs with showers and hadronization, leaving out the contributions 
from order a 2 s graphs, we shall refer to the program as LO+PS+Had. For comparison, we also 
display results from a pure lowest order perturbative calculation, LO, a pure next-to- leading 
order perturbative calculation, NLO, and just the Pythia program (version 6.221 with default 
parameters unless specially noted). Our standard choice for the cm. energy is ^/s = Mz 
with a s (M z ) = 0.118. In the perturbative calculations, we need to choose a renormalization 
scale. Our default choice is /i = y/s/6, based on the observation that in jet production for 
hadron physics the choice /z = Et/2 works well and in the e + e~ case each jet has a typical 
energy E ~ y/s/3. 

We will first examine the three-jet fraction / 3 as a function of the jet resolution pa- 
rameter y cut . Here the hope is that / 3 [NLO+PS+Had] will match / 3 [NL0] reasonably well 
since / 3 [NLO+PS+Had] is supposed to be correct to next-to- leading order. Then we will 
examine f^ 1 dfe/dM at a fixed choice for y cut , which we take as y cut = 0.05 (chosen be- 
cause at this value / 3 is approximately 1 x a s ). For / 3 _1 dfs/dM we recall that the pure 
NLO result is completely unrealistic. The hope is that f 3 _1 df 3 [NLO+PS+Had] /dM will match 
/3 -1 4^3 [Pythia] /dM reasonably well since Pythia is generally known to give a pretty realistic 
match to data. 

We begin with / 3 as a function of the jet resolution parameter y cut , which we display in 
Fig. 5. We first note that with our choice of renormalization scale, / 3 [L0] is rather close to 
/ 3 [NL0]. This agreement can be taken as a sign that the choice or renormalization scale was 
sensible. 

Next, we note that f 3 [LO+PS+Had] is quite a lot smaller than / 3 [NL0] or / 3 [L0]. What 
has happened is that the first stage of showering in the algorithm used here changes the 
normalization by an amount proportional to a s . This change is not as alarming as it might 
seem. For y cut = 0.05 we have / 3 [LO+PS+Had] // 3 [L0] 0.5. If we write this as one factor for 
each jet, 

fs [LO+PS+Had] r 

with a s (y/s/6) ~ 0.16, then the coefficient X is not large, X m 1.5. One could make X 
closer to zero by modifying the functions used in the primary showering, as discussed in the 
last paragraph of Sec. IV of Ref. [10]. For instance one could make the Sudakov exponent 
smaller. However, we do not attempt such an adjustment in this paper. 
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FIG. 5: Three-jet fraction, / 3 , versus y cut . The top panel shows / 3 [NL0]. The bottom panel 
shows the ratios of / 3 [NL0 + PS + Had], / 3 [L0 + PS + Had], / 3 [L0], and / 3 [Pythia] to / 3 [NL0]. The 
parameters are as in Figs. 1 and 2. 



We are now ready to look at / 3 [NL0+PS+Had], for which the correction terms based on 
order a 2 graphs are designed to produce a result that is correct to order a 2 . Indeed, we 
see that / 3 [NL0+PS+Had] matches / 3 [NL0] to about 10%, which is within the error that one 
expects for a next-to-leading order calculation. 

Finally in Fig. 5 we show / 3 [Pythia]. We see that / 3 [Pythia] matches / 3 [NL0] quite well, 
even though Pythia does not contain all the terms needed for next-to- leading order accuracy. 
On the other hand, Pythia does contain some next-to-leading order terms in the form of 
its choice of scale and in the generation of showers. Evidently, these terms do quite a good 
job of approximating the full next-to- leading order result. In fact, one might be inclined to 
rely entirely on Pythia if it were not for the fact that one cannot really tell how accurate 
its approximations are without having a full next-to-leading order calculation with which to 
compare it. 

A more rigorous perturbative test of whether NL0+PS, without hadronization in this case, 
is indeed correct to next-to-leading order can be obtained by studying the ratio 

= /3INL0+PS] JCjj 
' / 3 [»L0] ^ ' 

The difference between / 3 [NL0+PS] and / 3 [NL0] should be of order oP s so that the ratio R 
should have a perturbative expansion that begins at order a 2 . We can test this by evaluating 
/ 3 at different cm. energies yfs and plotting the ratio R against a 2 (y/s). The results of 
this test are shown in the upper panel of Fig. 6 for the cm. energies \fs = 34 GeV, Mz, 
500 GeV and 1000 GeV, corresponding to a s = 0.139, 0.118, 0.0940 and 0.0868, respectively. 
We see the expected shape of the R curve, which approaches a straight line through the 
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FIG. 6: Comparison of the NLO calculation with showers to a pure NLO calculation. In the top 
panel, we plot the ratio i2[NL0+PS] defined in Eq. (6). We take y cut = 0.05 and n = y/s/6. The ratio 
R is calculated for the cm. energies yfs = 34 GeV, Mz, 500 GeV and 1000 GeV, corresponding 
to a s = 0.139, 0.118, 0.0940 and 0.0868, respectively, and is plotted versus a 2 (^)/a 2 (M z ). For 
comparison, the curve R = 0.22a^(y / s) / 'ot 2 s {Mz) is shown. In the bottom panel, we show the 
difference between R calculated with hadronization and R calculated without hadronization. For 
comparison, we show the curve AR = 8 GeV/y/s. 



origin as a 2 s — > 0. For comparison, we show the line R = 0.22 a 2 {\fs) / a 2 {Mz)- In the 
lower panel of Fig. 6 we display the impact of hadronization by plotting the difference 
AR = i? [NLO+PS+Had] — i?[NL0+PS]. We observe that hadronization corrections reduce / 3 
by about 20% at ^Js = 34 GeV and by about 5% at -/s = Mz- They are negligible for 
larger cm. energies. For comparison, we show the curve AR = 8 GeV/yfs. The fact that 
the curve is similar to the data suggests that the combined NLO-MC program gives power 
law hadronization effects. 

A next-to-leading order calculation is supposed to do better than a leading order cal- 
culation because of its reduced uncertainty from uncalculated terms in the perturbative 
expansion. Some of the uncalculated terms contain logarithms of the renormalization scale. 
Thus one way to crudely estimate this uncertainty is to vary the renormalization scale by, 
say, a factor of two and see how the calculated quantity responds. In Fig. 7, we try this 
for / 3 [NL0+PS+Had] at y cnt = 0.05. We vary the renormalization scale \x in the perturba- 
tive part of the calculation while keeping the Pythia parameters unchanged. We see that 
/a [NLO+PS+Had] changes by about 5% when we increase or decrease \x by a factor 2. This is 
quite a small change that, we suspect, underestimates the uncertainty, since the difference 
between f 3 [NLO+PS+Had] and /afNLO] is 10% . The fj, dependence of the NLO calculation 
is even smaller. By way of comparison, we vary by a factor 2 the parameter in Pythia 
that controls the scale (PARJ(81)). As shown in Fig. 7, the result changes by about 15%, 
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consistent with the scale variation of the pure LO result. 
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FIG. 7: Three-jet fraction, /j versus n/no, where \x is the renormalization scale and (io is its default 
value. We show the result of the full NLO+PS+Had calculation and for the pure NL0 calculation. In 
these calculations the default scale is no = \/s/6. For comparison, we also show the result of 
varying the scale in plain Pythia. In this case, the default value, of the scale in a s ((i) is 
the transverse momentum in the parton splitting. The ratio fJ,/(J,o is the inverse of the Pythia 
parameter PARJ(81). We choose default parameters as in Fig. 1. 

Now we turn to the distribution of jet masses in three jet events, f^ 1 dfe/dM. We recall 
from Fig. 1 that f^ 1 df^\UUS\j dM is completely unrealistic in the small mass region. In 
Fig. 8, we plot dfe [NLO+PS+Had] jdM versus M. Both the unbounded increase in the 
cross section as M — > and the large negative contribution at M = are gone. In fact, 
the distribution looks a lot like the distribution produced by Pythia, Fig. 2, as can be seen 
from the direct comparison between / 3 -1 g?/ 3 [NLO+PS+Had] jdM and f^ 1 df 3 [Pythia] jdM in 
Fig. 9. It is worthy of notice in Fig. 9 that Pythia gets a result for / 3 _1 df^/dM that is in 
good agreement with df'3 [NLO+PS+Had] jdM without using next-to- leading order corrections. 

Finally, in Fig. 10 we investigate the effect of hadronization by compar- 
ing / 3 _1 df 3 [NLO+PS+Had] jdM to the same quantity with hadronization turned off, 
g(/3[NL0+PS] jdM. We see that without hadronization, both the unbounded increase in the 
NL0 cross section as M — > and the large negative contribution at M = are gone. How- 
ever, there is still a substantial probability to produce a jet with very small mass. With 
hadronization, it is no longer probable to produce a jet with mass less than 5 GeV. On the 
other hand, some 2 GeV is added to the mass of typical large mass jets, thus boosting the 
cross section at any fixed value of the jet mass above about 5 GeV. This effect explains most 
of the difference between / 3 _1 df 3 [NLO+PS+Had] / dM and rf/ 3 [NL0] jdM at large jet masses. 
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FIG. 8: Distribution of jet masses in three jet events, f^ 1 dfs/dM , in the full NLO+PS+Had calcu- 
lation using y cn t = 0.05. The result is substantially changed from Fig. 1 and is much closer to the 
result in Fig. 2. The calculation is defined as in Fig. 1. 
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FIG. 9: Distribution of jet masses in three jet events, f^ 1 dfs/dM , in the full NLO+PS+Had calcu- 
lation, in the Pythia calculation and in the pure NLO calculation. The calculations are defined as 
in Figs. 1 and 2. The NLO distribution is negative in the first mass bin, although this is not visible 
in a semilog plot. 
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FIG. 10: Comparison of / 3 _1 d/ 3 [NL0+PS+Had] jdM to the same quantity not including hadroniza- 
tion, / 3 _1 dfelNLO+PS] / dM . Hadronization increases the masses of jets. The calculation is defined 
as in Fig. 1. 

V. PARAMETER DEPENDENCE 

The algorithm for combining an NLO calculation with showers described in Refs. [10, 11] 
incorporates some adjustable parameters, most notably the two denoted as Ay and A so f t . 
The calculated value of an observable like the three-jet fraction is supposed to be correct 
to next-to-leading order independent of the choices of values of these parameters. Thus f% 
should be relatively insensitive Ay and A so ft- In this section we check whether this is so. 

We begin with Ay, which we can treat adequately with just two paragraphs. Consider 
the splitting of one of the partons from a Born graph, as represented by one of the square 
vertices in Fig. 3. We call this a primary splitting. In a primary splitting, there is an 
integration over the virtuality q 2 of the pair of daughter partons. In principle it is possible 
to use an integration range < q 2 < oo, since there is an automatic cutoff imposed by 
the kinematics of the calculation. However, one can impose a cut q 2 < Ay <f 2 , where q is 
the three-momentum carried by the pair of daughter partons. Taking Ay < oo changes the 
result from the Born graphs, but there is a compensating adjustment in the order af +l 
terms. One might imagine that taking Ay < oo is worthwhile because the parton splitting 
approximation with its Sudakov suppression, which is sensible at small virtualities, is then 
used only at small or moderate virtuality, while we use ordinary fixed order perturbation 
theory in the large virtuality region. The default value in the program is Ay = 1. We 
have calculated f% at y cut = 0.05 and yfs = M z for a range of values of Ay. We find that 
taking Ay < 0.5 makes / 3 [NL0+PS+Had] too big. However, in the range 0.5 < Ay < oo, f% is 
independent of Ay to within 10%. One might simply set Ay — > oo, thus eliminating it from 
the algorithm, except that very large values of Ay allow contributions from the region of 
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finite q 2 with q 2 — > 0. There is a singularity here (an artifact of the use of Coulomb gauge 
in the Feynman diagrams), producing undesirable fluctuations in the results. 

We examine next the parameter denoted as A so ft in Refs. [10, 11]. Consider the soft 
gluon emitted from the antenna produced by the three partons emitted in a Born graph, 
which is represented in Fig. 3 as a gluon emitted from the outgoing partons as a whole. 
The approximations used for the gluon emission are valid only when the gluon momentum 
is small. For this reason we impose a limit on the energy of the emitted gluon, 

E < M soft = A softv ^ (1 - *o) , (7) 

where y/s Q is the cm. energy of the three parton final state at the Born level (before show- 
ering) and t is the thrust value associated with this state. The default value of A so f t is 1/3. 
Changing A so ft changes the contribution from the Born graphs, but there is a compensating 
change in the order af +1 graphs, so that the net change is of order af +2 . 

In Fig. 11 we show how the calculated three-jet fraction depends on A so ft- We see that 



0.16 



0.15 



0.14 - 



0.13 



0.12 



e + e — > 3 jets: 3 -jet fraction f 3 

(Vs = M z ; li = Vs/6, k T algorithm, y = 0.05) 



NLO+PS+Had 



I 1 i 



0.2 



0.4 0.6 

^soft 



0.8 



FIG. 11: Dependence of [NLO+PS+Had] on the soft gluon cutoff parameter A so f t . The calculation 
is defined as in Fig. 1. 



varying A so ft over a substantial range, 0.1 < A so ft < 1, results in only a 10% change in the 
answer. We also see that, if we make A so ft much smaller than 1/3, the statistical uncertainty 
in the numerical integration increases. The reason is that in the af +1 graphs there are 
subtraction terms that remove the leading soft gluon singularities and the integration range 
covered by these subtraction is \k\ < M so ft. Effectively this inserts a cutoff \k\ > M so ft in the 
integrations over gluon momenta in the af +1 graphs. With a very small value of A so ft, the 
cutoff is removed. Then real emission graphs cancel virtual loop graphs after integration, 
but we encounter large positive and negative values point by point in the integration. 

There is another approach that one might have taken. One might leave the subtraction 
terms in the order af +1 graphs but simply omit the soft gluon emission from the Born graphs. 
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This is not valid at next-to-leading order when A so ft is finite, but it is allowed if A so ft is small 
enough. The reason is that the soft gluon emission from the Born graphs does not affect the 
measurement function as long as the gluon momentum is very small and the measurement 
function is infrared safe. Thus leaving out the soft gluon emission from the Born graphs is 
allowed if we check the dependence of the calculated quantity on A so ft in order to see how 
small is "small enough." This is similar to the procedure adopted in Refs. [6-8], where the 
corresponding cutoff parameter is called (3. We have tried this for the three-jet fraction and 
show the results in Fig. 12. What we find is that even A so ft = 1 is "small enough." That 
is, we need the subtraction terms, or some other soft gluon cutoff, in the af +1 graphs in 
order for the numerical integrations to work, but the three-jet fraction is so insensitive to 
soft gluons that one could do without the soft gluon emission that compensates for the net 
effect of the subtraction terms. This finding helps to justify the approach in Refs. [6-8]. 
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FIG. 12: Dependence of f% on the soft gluon cutoff parameter A so ft if the soft gluon emission from 
Born graphs is omitted. Also shown is f% [NLO+PS+Had] from Fig. 11. The calculation is defined as 
in Fig. 1. 



VI. CONCLUSIONS 

We have presented results from a method for adding parton showers to next-to-leading 
order calculations in QCD when the Born process involves massless strongly interacting 
partons. Specifically, the process investigated is e + e~ — > 3 jets. We start with an algorithm 
[10, 11] for coupling the NLO calculation to the first, hardest step in showering. Here the 
problem is to include the parton splittings that are part of the NLO calculation and the 
parton splittings that are part of the first step of showering without double counting. For 
our numerical results we have matched the NLO calculation to the standard Monte Carlo 
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event generator Pythia, although the method is more general and there is no limitation in 
principle to using another event generator. We have arranged for Pythia to be able to accept 
the partly developed shower and perform the rest of the showering plus hadronization. With 
an NLO program working together with Pythia, one can accomplish two things at once: 
calculate infrared-safe three jet quantities like the three jet fraction / 3 with the reduced 
theoretical error characteristic of a NLO calculation and get a sensible description of the 
inner structure of the three jets. 

Future NLO-MC hybrid programs will likely be more sophisticated than this one and 
less tied to a particular Monte Carlo program than that of Ref. [6-8]. They may employ 
the refined and simplified version of the matching algorithm presented in Ref. [16]. This 
newer algorithm is based on a commonly used subtraction formalism for performing NLO 
calculations, and it provides more flexibility by allowing to switch from a three-jet description 
to a two-jet description for events that are close to a two-jet configuration. 
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